use strict;
# 下载需要的文件

#&Download('https://www.ncbi.nlm.nih.gov/medgen/44376'); # MedGen上Metabolic disease的词条
#&Download('https://www.ncbi.nlm.nih.gov/medgen/6323'); # Inborn error of metabolism
&Download('https://ftp.ncbi.nih.gov/pub/medgen/MGSAT.RRF.gz'); # 这张表当中有Mesh转CUI的关系
&Download('ftp://nlmpubs.nlm.nih.gov/online/mesh/MESH_FILES/meshtrees/mtrees2018.bin');
&Download('ftp://nlmpubs.nlm.nih.gov/online/mesh/MESH_FILES/xmlmesh/desc2018.gz'); #这张表有疾病的详细信息
#&Download('https://ftp.ncbi.nih.gov/pub/medgen/medgen_pubmed_lnk.txt.gz');
# medgen到pubmed的对应关系。接近200M，但可从中提取MedGenID到Concept unique identifier（CUI），这个ID可以直接从Clinvar里的位点中筛选IMD相关的内容
&Download('https://ftp.ncbi.nih.gov/pub/clinvar/gene_condition_source_id'); # 基因与疾病的对应关系
# 7月5日的更新一定是临时工干的！目前clinvar的ftp上这个文件就剩下表头了
# 还好7月7号又重新更新了一下

&Download('https://ftp.ncbi.nih.gov/pub/clinvar/vcf_GRCh37/clinvar.vcf.gz'); # clinvar的vcf全部位点列表
&Download('https://ftp.ncbi.nih.gov/pub/clinvar/tab_delimited/allele_gene.txt.gz'); # 位点与基因的对应关系。既然clinvar有现成的，我不用自己搞这个了
ftp://nlmpubs.nlm.nih.gov/online/mesh/MESH_FILES/meshtrees/mtrees2018.bin
#system('gunzip -c desc2018.gz >desc2018.xml');
#system('Rscript MeshInfo.R'); # 利用R脚本，从xml文件当中抽出关键信息：第一列Mesh主ID，第二列为对应的概念名称，第三列为概念的层级信息

# #my $CMD='head -n 161 44376|tail -n 1|sed "s/></>\n</g"|grep  "/medgen/"';
# my $CMD='head -n 161 6323|tail -n 1|sed "s/></>\n</g"|grep  "/medgen/"';
# my $TXT=`$CMD`;
# my @a=split/\n/,$TXT;
# my %MedGenIDs=();
# my %IMDmedgenIDs=();
# my %IMDname=();
# for(my$i=1;$i<scalar(@a);$i++){
# 	# 跳过前5行。因为前面是Metabolic disease所属的大类，后边的才是各种属于IMD的疾病
# 	$a[$i]=~m/\/medgen\/(\d+)/;
# 	$MedGenIDs{$1}=1; # 记录IMD所属的MedGenID
# }
#
# open I,"gunzip -c medgen_pubmed_lnk.txt.gz|";
# while(<I>){
# 	chomp;
# 	my @a=split/\|/;
# 	if($MedGenIDs{$a[0]}){
# 		$IMDmedgenIDs{$a[1]}=$a[0];
# 		#print "$a[1]\n";
# 		$IMDname{$a[1]}=$a[2]; # 提取IMD的CUI，且提取IMD的英文名和medgen ID
# 	}
# }
# close I;

my $IMDmark='C16.320.565';
my %RelationInfo=();
my %IMDmesh=();
open I,"Mesh.xls" or die();
while (<I>) {
	unless(m/$IMDmark[\.\d+]+/){next;}
	@_=(m/($IMDmark[\.\d+]+)/);
	my @a=split/\t/;
	foreach my $x(@_){
		$RelationInfo{$x}=$a[0];
	}
	$IMDmesh{$a[0]}=$a[1];
}
close I;

my %IMDname=();
open I,"gunzip -c MGSAT.RRF.gz|";
while (<I>) {
	chomp;
	@_=split/\|/;
	unless($IMDmesh{$_[3]}){next;}
	$IMDname{$_[0]}=$_[3];
}
close I;

my %Selected=();
my @SelectedIDs=();
my %VarInfo=();
my %DiseaseAllele=();
my @Items=();
my $Head=(); # vcf文件表头
open I,"gunzip -c clinvar.vcf.gz|";
open O,"|bgzip -c >IMD.clinvar.vcf.gz"; # IMD专属VCF文件
while(<I>){
	if(m/\#\#/){
		print O $_;
		if(m/\#\#INFO\=\<ID\=(\S+?),/){
			if($1 eq 'CLNDISDB'){
				push(@Items,'MedGen');
				push(@Items,'MedGen_DiseaseName');
				push(@Items,'OMIM');
				push(@Items,'Orphanet');
				next;
			}
			if($1 eq 'CLNDN'){next;}
			if($1 eq 'CLNDNINCL'){next;}
			push(@Items,$1);
		}
		next;
	}
	if(m/^\#/){
		print O $_;
		chomp;
		s/\#//g;
		my @a=split/\t/;
		$Head=join("\t",@a[0..scalar(@a)-2]);
		next; # 提取总表的表头
	}
	chomp;


	my @a=split/\t/;
	my $clinvarID=$a[2]; # 这个是clinvar的主ID
	#print STDERR "$clinvarID\n";
	my @b=split/\;/,$a[-1];
	my $pass=0;
	foreach my $x(@b){
		my @c=split/\=/,$x;
		if($c[0] eq 'CLNDN'){next;} # 疾病名称我来统一撰写，不必用Clinvar自带的
		if($c[0] eq 'CLNDISDB'){
			my @u=split/\|/,$c[1]; #这么做是考虑单一位点关联多种疾病的情况
			foreach my $y(@u){
				my @v=split/\,/,$y;
				foreach my $z(@v){
					my @w=split/\:/,$z;
					$VarInfo{$w[0]}{$clinvarID}.="$w[1]|"; # 将多个数据库当中的疾病名称ID进行拆分存储。
					if(($w[0] eq 'MedGen')and($IMDname{$w[1]})){
						#print "$w[1]\n";
						$pass=1;
					} # 只要IMD的疾病ID被提取出来，则该位点就被选中了
				}
			}
			next;
		}
		$VarInfo{$c[0]}{$clinvarID}=$c[1];
	}
	if($pass){
		#print "$clinvarID\n";
		$Selected{$clinvarID}=join("\t",@a[0..scalar(@a)-2]); # 将选中的位点信息存起来
		push(@SelectedIDs,$clinvarID);
		print O "$_\n";
	}
}
close O;
close I;

my %GeneID=();
my %GeneSymbol=();
my %GeneName=();
open I,"gunzip -c allele_gene.txt.gz|";
open O,">IMD.allele_gene.xls"; # 位点到基因的映射文件，一度怀疑这张表真的是不是NCBI的临时工搞的……
$_=<I>;
print O $_;
while(<I>){
	my $clinvarID=(split/\t/,$_)[0];
	my @a=split/\t/;
	$GeneName{$a[1]}="$a[2]\t$a[3]"; # 基因名称
	unless($Selected{$clinvarID}){next;}
	unless($a[1] >0){next;} # 这份表格当中居然还有标记为-1的GeneID，舍掉它们
	$GeneID{$a[0]}.="$a[1]|";
	$GeneSymbol{$a[0]}.="$a[2]|";

	print O $_;
}
close I; # 打印一下致病位点所在的基因信息。可能有部分位点同时落在了多个基因上

# 接下来制作位点总表。首先包括VCF当中的关键内容，然后添加上位点所在的基因信息，尤其核查某些位点是不是没有落在基因上？
# my @Items=keys(%VarInfo); # 先把表头信息给我单独整出来
#print $Head;
my %Gene_Allele=();
my %Disease_Allele=();
$Head.="\t".join("\t",@Items);
open O,">IMD.allele.summary.xls";
print O "$Head\tGeneID\tGeneSymbol\n";
foreach (@SelectedIDs){
	print O $Selected{$_};
	foreach my $x(@Items){
		if($x eq 'MedGen'){
			my $txt='-';
			@_=();
			if($VarInfo{'MedGen'}{$_}){
				foreach my $y(split/\|/,$VarInfo{'MedGen'}{$_}){
					unless($IMDname{$y}){next;}
					$Disease_Allele{$y}++; #统计每种疾病中突变位点的数目
					push(@_,$y);
				}
				$txt=join("|",@_);
			}
			print O "\t$txt";
			next;
		}
		if($x eq 'MedGen_DiseaseName'){
			my $txt='-';
			@_=();
			if($VarInfo{'MedGen'}{$_}){
				foreach my $y(split/\|/,$VarInfo{'MedGen'}{$_}){
					unless($IMDname{$y}){next;}
					push(@_,$IMDname{$y});
				}
				$txt= join("|",@_);
			}
			print O "\t$txt";
			next;
		}
		if($x eq 'RS'){
			my $txt='-'; # 空白占位符
			if($VarInfo{'RS'}{$_}>0){
				$txt="rs$VarInfo{RS}{$_}";
			}

			print O "\t$txt";
			next;
		}
		my $txt='-'; # 空白占位符
		if($VarInfo{$x}{$_}){
			@_=split/\|/,$VarInfo{$x}{$_};
			$txt=join("|",@_);
		}
		print O "\t$txt";
	}

	if($GeneID{$_}){
		@_=split/\|/,$GeneID{$_};
		foreach my $z(@_){ # 数数一个基因上多少位点
			$Gene_Allele{$z}++;
		}
		my $txt= join("|",@_);
		print O "\t$txt";
	}else{
		print O "\t-";
	}

	if($GeneSymbol{$_}){
		@_=split/\|/,$GeneSymbol{$_};
		my $txt= join("|",@_);
		print O "\t$txt";
	}else{
		print O "\t-";
	}

	print O "\n";
}
close O;
close I;

my %Gene2Disease=();
open I,"<gene_condition_source_id";
open O,">IMD.gene_condition.xls";
$Head=<I>;
print O $Head;
while(<I>){
	my @a=split/\t/;
	unless($IMDname{$a[3]}){next;}
	#print "$a[3]\n";
	print O $_;
	$Gene2Disease{"$a[0]\t$a[3]"}=1; # 记录基因与疾病的对应关系，待会儿来解析
}
close O;
close I;

open O,">IMD.gene.summary.xls"; # IMD相关基因的总表
print O "GeneID\tGeneSymbol\tGeneName\tMedGen_CUI\tMedGen_Disease\tDisease_Count\tAllele_Count\n";
my %temp=();
foreach (keys %Gene2Disease){
	my @a=split/\t/;
	$temp{$a[0]}.="$a[1]|";
}
foreach (keys %temp){
	my @a=split/\|/,$temp{$_};
	my @b=@a;
	my $Count = scalar(@a);
	foreach (my $i=0;$i<scalar(@a);$i++){
		$b[$i]=$IMDname{$a[$i]};
	}
	my $MGid=join("|",@a);
	my $MGname=join("|",@b);
	unless($Gene_Allele{$_}){$Gene_Allele{$_}=0;}
	print O "$_\t$GeneName{$_}\t$MGid\t$MGname\t$Count\t$Gene_Allele{$_}\n";
}
close O; # 经过检查，发现clinvar给出的部分与疾病相关的基因中，并没给出基因上的突变位点。推测可能是突变位点调控了某基因，但却不在哪个基因区间当中

open O,">IMD.disease.summary.xls"; # IMD疾病总表
print O "MedGen_CUI\tDiseaseName\tMedGen_Disease\tGene\tGene_Count\tAllele_Count\n";
%temp=();
foreach(keys %Gene2Disease){
	my @a=split/\t/;
	$temp{$a[1]}.="$a[0]|";
}
foreach (keys %temp){
	my @a=split/\|/,$temp{$_};
	my @b=@a;
	my $Count=scalar(@a);
	foreach(my $i=0;$i<$Count;$i++){
		$b[$i]=(split/\t/,$GeneName{$a[$i]})[0];
	}
	my $GNames=join("|",@b);
	print O "$_\t$IMDmesh{$IMDname{$_}}\t$IMDname{$_}\t$GNames\t$Count\t$Disease_Allele{$_}\n";
}
close O;

sub Download{
  # 若文件存在，就不必下载了。若文件不在则自动下载文件
  my $link=shift;
  my $filename=(split/\//,$link)[-1];
  unless(-s $filename){
    system("wget $link");
  }
}
